import numpy as np
from sympy import Matrix
from scipy import linalg
from fractions import Fraction
import sympy as sp
import math
from math import exp
#1.求Jordan标准性-start
a = np.array([[4,3,0,1], [0,2,0,0], [1,3,2,1],[0,0,0,2]])#-------------输入-------------
# a = np.array([[2*3.14,0,0,0,0,0],[0,2*3.14,2,0,0,0],[0,0,2*3.14,0,0,0],
#               [0,0,0,-2*3.14,2,0],[0,0,0,0,2*3.14,2],[0,0,0,0,0,-2*3.14,]])
# a = np.array([[3,-7,-6],
#               [2,-2,-2],
#              [-2,2,2]])
# print(-1*a.dot(a))
print("a:\n",a)
m = Matrix(a)
P, J = m.jordan_form()#P-相似变换矩阵，J-Jordan标准形（可看出特征值及其几何重数）
P=np.array(P)
P_inv=linalg.pinv(P.astype(np.float64) )#P的逆矩阵-先将MutableDenseMatrix转换为array
np.set_printoptions(precision=2) #设置np.array的小数位数
# np.set_printoptions(formatter={'all':lambda x: str(Fraction(x).limit_denominator())})#设置np.array用分数表示
print("J-Jordan标准形:",J)
#1.求Jordan标准性-end
print("相似变换矩阵P:",P)
print("逆矩阵P-1:",P_inv)#P的逆矩阵
# rref = Matrix(np.array(P)).rref()[0].tolist()#矩阵P的行最简式
# print("矩阵P的行最简式：",rref)

print("矩阵的秩:",np.linalg.matrix_rank(a))#返回矩阵的秩
np.linalg.det(a) #返回矩阵的行列式
a.diagonal() #返回矩阵的对角线元素，也可以通过offset参数在主角线的上下偏移，获取偏移后的对角线元素。a.diagonal(offset=1)返回array([1.10])
a.trace()#返回迹,主对角线上各个元素的总和被称为矩阵A的迹（或迹数），一般记作tr(A)。
eigenvalues ,eigenvectors= np.linalg.eig(a) #eigenvalues 为特征值。eigenvectors为特征向量#要求是方阵
print("eigenvalues:",eigenvalues)
for λ in eigenvalues:
    eigen_matrix = a - np.eye(a.shape[0]) * λ  # 特征矩阵A-λE(对指定特征值)
    #控制输出：list保留小数位数
    print("特征值", '%.2f'%λ, "有", a.shape[0] - np.linalg.matrix_rank(eigen_matrix), "个线性无关的特征向量")  # 特征矩阵的秩
# λ=int(eigenvalues[1])#取特征值
# eigen_matrix=a-np.eye(a.shape[0])*λ#特征矩阵A-λE(对指定特征值)
print("eigenvectors:",eigenvectors)

#-------------计算m阶Jordan函数f(J)和方阵函数f(A)-------------------
x,y,t= sp.symbols('x y t')
# func = sp.exp(x*t)##-------------输入-------------幂函数,  # 多变量符号函数(带变量符号)
func = sp.exp(x)##-------------输入-------------幂函数,f(x)=e^x  #符号函数(带入值)
# func = sp.sin(x)#正弦函数,f(x)=sin(x)
# func = sp.cos(x)#余弦函数,f(x)=cos(x)
x1 = eigenvalues[0]  # x-取特征值

eig_set=set(eigenvalues)#特征值集合
print("特征值=",x1,",函数f(x)=",func)#

f_J_diag_temp=[]
for eig in eig_set:
    f_J_λ = []  # m阶Jordan块矩阵-某个特征值的
    # print(i)
    r = list(eigenvalues).count(eig)  # 计算特征值阶数
    for i in range(r):
        tmp = []
        for j in range(r):
            if j<i:
                tmp.append(0)
            else:
                tmp_der=sp.diff(func,x,j-i)#符号函数求导
                # func_der1_x1=tmp_der/(math.factorial(j-i))# # 符号函数(带变量符号)
                func_der1_x1 = float(tmp_der.evalf(subs={x: eig}))/(math.factorial(j-i))  # 符号函数(带入值)-用值替换符号变量x
                tmp.append(func_der1_x1)
        f_J_λ.append(tmp)
    f_J_diag_temp.append(f_J_λ)
    print("特征值=",eig,"----f_J_λ:",f_J_λ)
#由list构建对角块矩阵
f_J_diag=[]
for k in range(len(f_J_diag_temp)):
    f_J_diag = linalg.block_diag(f_J_diag,f_J_diag_temp[k])
f_J_diag=np.array(f_J_diag[1:])#得到m阶Jordan块矩阵-首行为空，切片掉
print("m阶Jordan块矩阵f_J_diag:\n",f_J_diag)
print("f(x)=",func)
f_A_=np.dot(np.dot(P,f_J_diag),P_inv)#计算方阵函数"f(A)=P*f(J)*P-1:"
print("方阵函数-带入值f(A)=P*f(J)*P-1:\n",np.array(f_A_))
print("---------------------------------")
#-------------计算m阶Jordan函数f(J)和方阵函数f(A)-------------------

#-------------计算m阶Jordan函数f(J)和方阵函数f(A),带变量符号-------------------
x,y,t= sp.symbols('x y t')
func = sp.exp(x*t)##-------------输入-------------幂函数,  # 多变量符号函数(带变量符号)
# func = sp.exp(x)##-------------输入-------------幂函数,f(x)=e^x  #符号函数(带入值)
# func = sp.sin(x)#正弦函数,f(x)=sin(x)
# func = sp.cos(x)#余弦函数,f(x)=cos(x)
x1 = eigenvalues[0]  # x-取特征值

eig_set=set(eigenvalues)#特征值集合
print("特征值=",x1,",函数f(x)=",func)#

f_J_diag_temp=[]
for eig in eig_set:
    f_J_λ = []  # m阶Jordan块矩阵-某个特征值的
    # print(i)
    r = list(eigenvalues).count(eig)  # 计算特征值阶数
    for i in range(r):
        tmp = []
        for j in range(r):
            if j<i:
                tmp.append(0)
            else:
                tmp_der=sp.diff(func,x,j-i)#符号函数求导
                func_der1_x1=tmp_der.subs(x,eig)/(math.factorial(j-i))# # 符号函数(带变量符号)
                # func_der1_x1 = float(tmp_der.evalf(subs={x: eig}))/(math.factorial(j-i))  # 符号函数(带入值)-用值替换符号变量x
                tmp.append(func_der1_x1)
        f_J_λ.append(tmp)
    f_J_diag_temp.append(f_J_λ)
    print("特征值=,",eig,"----f_J_λ:",f_J_λ)
#由list构建对角块矩阵
f_J_diag=[]
for k in range(len(f_J_diag_temp)):
    f_J_diag = linalg.block_diag(f_J_diag,f_J_diag_temp[k])
f_J_diag=np.array(f_J_diag[1:])#得到m阶Jordan块矩阵-首行为空，切片掉
print("m阶Jordan块矩阵f_J_diag:\n",f_J_diag)
print("f(x)=",func)
f_A_=np.dot(np.dot(P,f_J_diag),P_inv)#计算方阵函数"f(A)=P*f(J)*P-1:"
print("方阵函数-带变量符号,f(A)=P*f(J)*P-1:\n",np.array(f_A_))
#-------------计算m阶Jordan函数f(J)和方阵函数f(A),带变量符号-------------------